Exploring New Alloy Systems with Pymatgen¶
Outline¶
-
Select a test-case system
- 1.1 Exercise:
StructureandMPResterrefresher - 1.2 Lesson: add oxidation states to a
Structure - 1.3 Bonus: plot
BandStructure
- 1.1 Exercise:
-
Select an alloy partner
- 2.1 Lesson: find possible dopants
- 2.2 Exercise: find the best alloy partner (A = ?) for AxZn1-xS
- 2.3 Lesson: explore phase diagrams
-
Transform to make a new CuxZn1-xS alloy
- 3.1 Lesson: structure transformation
- 3.2 Exercise: try your own transformation on CuZnS2
-
Calculate new properties
- 4.1 Lesson: volume prediction and XRD plot
- 4.2 Exercise: try this on your CuZnS2 structure
-
Test your skills
- 5.1 Exercise: compare relaxed DFT structures to estimates
- 5.2 Lesson: add computed entries to phase diagram
- 5.3 Next steps
1. Select a test-case system¶
In this notebook we will focus on cubic zinc-blende ZnS, a wide band gap (transparent) semiconductor. In my PhD research I study p-type transparent semiconductors, so I will pose the question: how can we use ZnS as a starting point to create a p-type transparent semiconductor, and how can pymatgen help with this?
Import the MPRester client:
from pymatgen.ext.matproj import MPRester
The Materials ID (mp-id) of zinc-blende ZnS is mp-10695, see https://materialsproject.org/materials/mp-10695/.
ZnS_mpid = "mp-10695"
1.1 Exercise: Structure and MPRester refresher¶
Get the structure¶
with MPRester() as mpr: # YOUR API KEY GOES IN THIS FUNCTION! copy from materialsproject.org/dashboard
ZnS_structure = mpr.get_structure_by_material_id(ZnS_mpid)
#### if you're having problems with your internet or API key:
# from monty.serialization import loadfn
# ZnS_structure = loadfn("assets/ZnS_structure.json")
Get space group information¶
ZnS_structure.get_space_group_info()
If you want to, try it out on our web app here.
- Click "Draw atoms outside unit cell bonded to atoms within unit cell"
- Play around with it!
1.2 Lesson: add oxidation states to a Structure¶
Pymatgen has a simple transformation to estimate the likely oxidation state of each specie in stoichiometric compounds using a bond-valence analysis approach. This information is needed to compare ionic radii and assess substitutional dopant probability. You can also enter the oxidation states manually if you'd prefer.
from pymatgen.transformations.standard_transformations import AutoOxiStateDecorationTransformation
Initialize this transformation:
oxi_transformation = AutoOxiStateDecorationTransformation()
ZnS_structure_oxi = oxi_transformation.apply_transformation(ZnS_structure)
print(ZnS_structure_oxi)
1.3 Bonus: plot BandStructure¶
from pymatgen.electronic_structure.plotter import BSPlotter
This code retrieves a BandStructureSymmLine object which contains all the information about a line-mode band structure.
with MPRester() as mpr: # YOUR API KEY GOES IN THIS FUNCTION! copy from materialsproject.org/dashboard
ZnS_bs = mpr.get_bandstructure_by_material_id(ZnS_mpid)
#### if you're having problems with your internet or API key
# from monty.serialization import loadfn
# ZnS_bs = loadfn("assets/ZnS_bs.json")
This band structure can be plotted using BSPlotter:
ZnS_bsp = BSPlotter(ZnS_bs)
ZnS_bsp.show() # takes a second
Band gap correction¶
ZnS_bs.get_band_gap()
ZnS has an experimental gap of approximately 3.5 eV, but the GGA calculated gap is far too low! We can apply a "scissor" to this band structure to correct for this.
Scissor corrections are only appropriate if they're clearly acknowledged! They're used here because we expect the shape of the bands to be correct from our GGA calculation, and we know experimentally the gap is 3.5 eV. We do not use any scissor corrections on the Materials Project website or in the database.
ZnS_bs_scissor = ZnS_bs.apply_scissor(new_band_gap=3.5)
ZnS_bsp_scissor = BSPlotter(ZnS_bs_scissor)
ZnS_bsp_scissor.show()
2. Select an alloy partner¶
2.1 Lesson: find possible dopants¶
Scientific question: Which p-type dopants are most likely to sit at substitutional sites in ZnS?
Pymatgen has a machine-learned method for estimating the probability that one ion will substitute for another (Hautier et al. 2011), and reports the results ranked in order of probability. Note the input structure has to be "decorated" with oxidation states for this method to work.
from pymatgen.analysis.structure_prediction.dopant_predictor import get_dopants_from_substitution_probabilities
substitutional_dopants = get_dopants_from_substitution_probabilities(
ZnS_structure_oxi, num_dopants=10)
Here are some options to dope ZnS p-type:
p_dopants = substitutional_dopants['p_type']
We can see this returns a list of dictionaries:
print(p_dopants)
To make this easier to read we can use the pandas package:
import pandas as pd
pd.DataFrame(p_dopants)
2.2 Exercise: find the best alloy partner (A = ?) for AxZn1-xS¶
Scientific question: is a p-type zinc-blende AxZn1-xS alloy possible?
Let's see if zinc-blende binaries exist for these ternaries, and how far off the hull they sit.
Find dopants¶
First, find a list of possible cation dopant elements:
# I've pre-written this code block for convenience
# all it does is take the possible dopants list given previously, takes the cations, and makes a list of their elements
possible_cation_dopants = []
for x in p_dopants:
specie = x["dopant_species"]
if specie.oxi_state > 0:
possible_cation_dopants.append(str(specie.element))
print(possible_cation_dopants)
Query for end-point structure¶
Next, let's query the MPRester to make a table of all of the binary compounds with a space group "F-43m" that contain sulfur and one of these possible_cation_dopants. Note that the query criteria are listed on the mapidoc.
with MPRester() as mpr: # YOUR API KEY GOES IN THIS FUNCTION! copy from materialsproject.org/dashboard
query = mpr.query(
{
# the query criteria
"elements": {
"$all": ["S"], "$in": possible_cation_dopants
},
"nelements": 2,
"spacegroup.symbol": "F-43m"
},
# the properties we want to return
[
"task_id",
"e_above_hull",
"pretty_formula",
"theoretical",
"spacegroup.symbol"
]
)
#### if you're having problems with your internet or API key
# from monty.serialization import loadfn
# query = loadfn("assets/alloy_partner_query.json")
pd.DataFrame(query)
Which cation should we pick? Cu! Ag8S a theroetical intermetallic, and it's energy is ridiculously high. CuS is a theoretical compound and is not "on the hull," but it's close at only 0.01 eV/atom, meaning it is only slightly metastable. Ok, so let's pick Cu+ to use as a p-type dopant, which I've explored experimentally in the past (see Woods-Robinson et al. 2019).
Retrieve end-point structure¶
To proceed, we have to retrieve the Structure for CuS:
CuS_mpid = query[0]["task_id"]
with MPRester() as mpr: # YOUR API KEY GOES IN THIS FUNCTION! copy from materialsproject.org/dashboard
CuS_structure = mpr.get_structure_by_material_id(CuS_mpid)
#### if you're having problems with your internet or API key
# from monty.serialization import loadfn
# CuS_structure = loadfn("assets/CuS_structure.json")
Yep! We’re not done, but this is a good starting point for dopants to investigate with further defect calculations. This can be accomplished using workflows from packages like PyCDT (Broberg et al. 2018) which integrate with pymatgen's defect capabilities.
2.3 Lesson: explore phase diagrams¶
Scientific question: what does Cu-Zn-S phase space look like?
There are many built-in tools to explore phase diagrams in pymatgen. To build a phase diagram, you must define a set of ComputedEntries with compositions, formation energies, corrections, and other calculation details.
from pymatgen.analysis.phase_diagram import PhaseDiagram, PDPlotter, CompoundPhaseDiagram, GrandPotentialPhaseDiagram
We can import entries in this system using the MPRester. This gives a list of all of the ComputedEntries on the database:
with MPRester() as mpr: # YOUR API KEY GOES IN THIS FUNCTION! copy from materialsproject.org/dashboard
entries = mpr.get_entries_in_chemsys(['Cu', 'Zn', 'S'])
#### if you're having problems with your internet or API key:
# from monty.serialization import loadfn
# entries = loadfn("assets/Cu-Zn-S_entries.json")
phase_diagram = PhaseDiagram(entries)
Conventional phase diagram¶
plotter = PDPlotter(phase_diagram, show_unstable=True, markersize=20) # we increase the marker size here to make it easier to see the stable points
plotter.show()
And check the space group:
CuZnS2_ordered_structure.get_space_group_info()
Is it the same as ZnS? Because of the Cu substitution, this structure has a different space group than ZnS!
4. Calculate new properties¶
4.1 Lesson: volume prediction and XRD plot¶
So far we just have a really rough guess of an alloy structure, and the lattice parameters are still equal to those of ZnS. We can estimate the new volume \(V_{x-guess}\) after the substitution using Vegard's Law (assuming zero bowing).
$V_{x-estimate} = V_{scaling } \times [ V_{CuS}(x) + V_{ZnS}(1-x) ] $
$V_{CuZn_3S_4-estimate} = [2\times2\times2] \times [ V_{CuS}(0.25) + V_{ZnS}(0.75) ] $
scaling_matrix
x
scaling_volume = scaling_matrix.prod()
CuZn3S4_estimated_volume = scaling_volume * ((ZnS_structure.volume) * (1 - x) +
(CuS_structure.volume) * x)
print(CuZn3S4_ordered_structure.volume)
print(CuZn3S4_estimated_volume)
CuZn3S4_structure_estimate = CuZn3S4_ordered_structure.copy()
CuZn3S4_structure_estimate.scale_lattice(CuZn3S4_estimated_volume)
print(CuZn3S4_structure_estimate)
This is better but still wrong, and does not take into account any structural distortions. Note that there are some other methods on pymatgen to guess structure volume (see pymatgen.analysis.structure_prediction.volume_predictor), but in my experience Vegard's law is usually just as helpful. Your next step would be to relax this new structure using DFT or another method (see below).
Calculate XRD, compare to original structure¶
Now we can compare this structure to our original ZnS and CuS structure to, for example, see how the X-ray diffraction (XRD) patterns are expected to shift as x increases in CuxZn1-xS:
from pymatgen.analysis.diffraction.xrd import XRDCalculator
Initialize the XRDCalculator with the conventional Cu-K\(\alpha\) wavelength (note: Cu here has nothing to do with the Cu we're adding to the structure):
xrd = XRDCalculator()
structures = [
ZnS_structure,
CuZn3S4_structure_estimate,
CuS_structure
]
xrd_plots = xrd.plot_structures(
structures,
two_theta_range=[25, 65],
annotate_peaks=False, # to keep the plot cleaner
size_kwargs={'w': 10, 'h': 9}, # these options are optional to make the plot look nicer
)
You can see how the \(2\theta\) peaks shift slightly to the right with addition of Cu!
4.2 Exercise: try this on your CuZnS2 structure¶
Guess the structure volume using Vegard's Law, and correct for this:
$V_{x-estimate} = V_{scaling } \times [ V_{CuS}(?) + V_{ZnS}(?) ] $
x_CuZnS2
scaling_volume
CuZnS2_structure_estimate = CuZnS2_ordered_structure.copy()
CuZnS2_structure_estimate.scale_lattice(scaling_volume *
((ZnS_structure.volume) * (1 - x_CuZnS2) +
(CuS_structure.volume) * x_CuZnS2))
Print the new structure
print(CuZnS2_structure_estimate)
Add this structure to the series of XRD plots to compare XRD for x = 0, 0.25, 0.5, 1:
structures = [
ZnS_structure,
CuZn3S4_structure_estimate,
CuZnS2_structure_estimate,
CuS_structure
]
xrd_plots = xrd.plot_structures(
structures,
annotate_peaks=False,
two_theta_range=[25,65],
size_kwargs={'w': 10, 'h': 12}
)
5. Test your skills¶
This is the wee beginning of making an alloy. Here are some follow-up steps:
I constructed similar alloys to those that we just explored, at x = 1/4 and x = 1/2, and relaxed them with DFT. We'll explore my results here:
5.1 Exercise: compare relaxed DFT structures to estimates¶
from pymatgen.core.structure import Structure
These are my output .cif files from one of our DFT workflows. See fireworks and atomate packages for details here.
CuZn3S4_relaxed = Structure.from_file("assets/Zn3CuS4_Amm2.cif")
CuZnS2_relaxed = Structure.from_file("assets/ZnCuS2_R3m.cif")
How do these space groups compare to our estimates?¶
CuZn3S4_structure_estimate.get_space_group_info()
CuZn3S4_relaxed.get_space_group_info()
CuZnS2_structure_estimate.get_space_group_info()
CuZnS2_relaxed.get_space_group_info()
Are they higher or lower in symmetry? After relaxation, both structures are in a different space group (lower symmetry) than the alloys we just made. This is likely due to structural distortions.
Add in DFT structure to XRD¶
Replace the alloy structures in the previous XRD exercise with the two relaxed alloy structures, again comparing XRD for x = 0, 0.25, 0.5, 1:
structures = [
ZnS_structure,
CuZn3S4_relaxed,
CuZnS2_relaxed,
CuS_structure
]
xrd_plots = xrd.plot_structures(
structures,
annotate_peaks=False,
two_theta_range=[25,65],
size_kwargs={'w': 10, 'h': 12},
)
Peak splittings are now present in the diffraction patterns, and the shift to higher \(2\theta\) is not as significant.
5.2 Lesson: add computed entries to phase diagram¶
Scientific question: are these new phases stable?
To assess the stability of these new phases, let's look at JSON files containing ComputedEntry data:
from monty.serialization import loadfn
Zn3CuS4_Amm2_entry = loadfn("assets/Zn3CuS4_Amm2_entry.json")
ZnCuS2_R3m_entry = loadfn("assets/ZnCuS2_R3m_entry.json")
These entries were created by relaxing the above structure using one of our DFT workflows. An "entry" is mainly just a composition and an energy, so can be created manually, without performing a calculation, or even from experimental data
print(ZnCuS2_R3m_entry)
We can add these two entries to entries, our set of ComputedEntry data from MP in the Cu-Zn-S phase space:
new_entries = entries + [Zn3CuS4_Amm2_entry, ZnCuS2_R3m_entry]
new_phase_diagram = PhaseDiagram(new_entries)
Conventional phase diagram¶
new_plotter = PDPlotter(new_phase_diagram, show_unstable=10, markersize=20)
x = new_plotter.get_plot(new_phase_diagram, label_unstable=False)
We see our two new phases show up here! How does the energy landscape change?
Contour phase diagram¶
new_fig = new_plotter.get_contour_pd_plot()
Compare to the phase diagram before new phases were added:
fig = plotter.get_contour_pd_plot()
Binary phase diagram¶
from pymatgen.core.composition import Composition
new_cpd = CompoundPhaseDiagram(new_entries, [Composition("ZnS"), Composition("CuS")])
new_compound_plotter = PDPlotter(new_cpd, show_unstable=10, markersize=20)
new_compound_plotter.show()
new_phase_diagram.get_e_above_hull(ZnCuS2_R3m_entry)
new_phase_diagram.get_e_above_hull(Zn3CuS4_Amm2_entry)
Because these phases break the binary hull, this shows that there is likely a stable phase of CuZnS2, though this may not be the lowest energy phase. Zn3CuS4 is highly metastable.
5.3 Next steps¶
Here are some follow-up calculations you can do to explore ternary space:
-
Structure Prediction: To explore other possibilities of polymorphs in Cu-Zn-S ternary space, one could perform a "structure prediction" in this chemical space. You can use the Materials Project website's Structure Prediction app. There are also methods in
pymatgen.analysis.structure_predictionto do this. -
DFT Calculations: You can submit your structure(s) to MP by loading it with the Crystal Toolkit app and clicking "Submit to Materials Project". You can use
fireworksandatomateto submit DFT calculations yourself to relax the structure, refine it, calculate band structure, and other properties. -
Local Environment Analysis: With your relaxed structures, you can upload your structure to the Crystal Toolkit app and look at "Analyze: Local Environments" to compare coordination environments. You can also explore the features in
pymatgen.analysis.local_env. -
Substrate Matcher: If you want to grow this alloy structure epitaxially, you can explore different substrates to use on the website using methods in
pymatgen.analysis.substrate_analyzer module.